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Abstract. The effects of self-consistently including superthermal electrons in the definition of 
the ambipolar electric field are investigated for the case of plasmaspheric refilling after a geomag- 
netic storm. By using the total electron population in the hydrodynamic equations, a method for 
incorporating superthermal electron parameters in the electric field and electron temperature calcu- 
lation is developed. Also, the ambipolar electric field is included in the kinetic equation for the 
superthermal electrons through a change of variables using the total energy and the first adiabatic 
invariant. Calculations based on these changes are performed by coupling time-dependent models 
of the thermal plasma and superthermal electrons. Results from this treatment of the electric field 
and the self-consistent development of the solution are discussed in detail. Specifically, there is a 
decreased thermal electron density in the plasmasphere during the first few minutes of refilling, a 
slightly accelerated proton shock front, and a decreased superthermal electron flux due to the decel- 
eration by the electric field. The timescales of plasmaspheric refilling are discussed and determined 
to be somewhat shorter than previously calculated for the thermal plasma and superthermal 
electron population due to the effects of the field-aligned potential. 


1. Introduction 

The plasmasphere is the region of near-Earth space where 
geomagnetic flux tubes corotate with the Earth [Nishida, 
1966], The convection electric field sweeps away the flux 
tubes outside the plasmasphere, never allowing them to fill to 
a steady state level of thermal plasma. During a geomagnetic 
disturbance, the plasmasphere is reduced due to the increased 
convection, sweeping away the plasma in the flux tubes that 
were recently corotating. When the convection decreases, 
these flux tubes return to their corotating trajectories around 
the Earth, and refilling of the thermal plasma can begin (see 
reviews by Horwitz [1987]; Singh and Horwitz [1992]; and 
Singh et al [1994]). 

During refilling after a geomagnetic storm, the ions can 
flow from the ionosphere into the plasmasphere at supersonic 
speeds [Banks et al ., 1971; Khazanov et al. , 1984; Singh et 
al, 1986; Rasmussen and Sc hunk, 1988; Guiter and Gombosi , 
1990; Gorbachev et al. , 1991; Wilson et al., 1992; Guiter et 
al. y 1995], resulting in shock formations. The flow condi- 
tions vary dramatically along a flux tube, from a subsonic, 0 + 
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dominated, collisional region at low altitudes to a supersonic, 
H + dominated, collisionless region at high altitudes. In the 
high-altitude region, highly non-Maxwellian distributions 
develop as the particles stream into the plasmasphere. The 
high mobility of the electrons allows them to outrun the ions; 
this violates charge neutrality, however, and an electric field 
appears to force the electrons to move with the ions (see 
Figure 1). This drag on the electrons is also an acceleration 
mechanism for the ions, which changes the potential struc- 
ture. 

Superthermal electrons are created in the ionosphere by 
photoionization and impact ionization of atmospheric neu- 
trals. These electrons, which have a highly structured source 
function and represent a nonthermal tail in the electron distri- 
bution function, must be treated kinetically in plasma models 
[Khazanov et al, 1994]. They are affected by the electric field 
and the accelerated ions and can redistribute the potential 
along the field line, and a self-consistent calculation should be 
performed in order to model the dynamics of this phenomenon 
correctly. 

A recent approach to coupling the thermal plasma popu- 
lations is to calculate the electron density and velocity from 
the conditions of quasi-neutrality and current balance, and then 
find the parallel electric field and the electron temperature from 
the electron momentum and energy equations. These values 
are then used in the equations for the other species. This type 
of calculation was performed for the inner magnetosphere 
[Khazanov et al, 1984; Richards and Torr, 1986; Guiter and 
Gombosi, 1990; Gorbachev et al, 1991; Miller et al, 1993; 
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Figure 1 . Schematic of a plasmaspheric flux tube, showing 
the helical path of a superthermal electron and the usual direc- 
tion of the ambipolar electric field. 


Singh et al., 1994; Gutter et al. 1995], the topside ionosphere 
in the auroral region [Min et al., 1993], and the polar iono- 
sphere [Schunk and Watkins , 1981; Gombosi and Nagy, 1989; 
Tam et al., 1995]. 

Several of these models calculate superthermal electrons and 
the thermal plasma simultaneously at various levels of self- 
consistency. The model of Khazanov et al. [1984] calculates 
the thermal plasma with a time-dependent 16-moment trans- 
port model for the electrons and 0 + and H + ions, along with a 
steady state, superthermal electron kinetic model for thermal 
heating rates. The field line interhemispheric plasma (FLIP) 
model [Richards and Torr , 1986] calculates the thermal plasma 
with a time-dependent five-moment fluid model and uses a 
steady state two-stream model for the superthermal electrons 
[e.g., Nagy and Banks , 1970] for the heating rates. The model 
of Gorbachev et al. [1991] includes wave activity in its 13- 
moment equation set for the thermal plasma calculation, 
obtaining slightly different results, especially in the morn- 
ing/evening sectors when refilling/depletion of the flux tube 
occurs. Min et al. [1993] calculate E\\ from the thermal plasma 
equations and then use this electric field in the steady state 
kinetic equation for the superthermal electrons in the aurora 
and at midlatitudes. Most recently, Tam et al. [1995] treat the 
ions and superthermal electrons kinetically along with a fluid 
approach for the thermal electrons, and obtain steady state 
polar wind results that are collisionally and electrodynam- 
ically self-consistent. 

The goal of this study is to extend results such as these to a 
time-dependent calculation while collisionally and electrody- 
namically coupling the plasma populations. This will 
involve a reexamination of the coupling processes between 
the superthermal electrons and the thermal plasma, resulting 
in various terms linking the equation sets. A driver program is 
used that can couple any thermal plasma model to the 
superthermal electron model of Khazanov and Liemohn 
[1995]. For this study, the two-stream, field-aligned, five- 
moment hydrodynamic model of Guiter et al. [1995] will be 
used for the thermal plasma calculation. Since plasmaspheric 
refilling along closed field lines is the focus of this study, this 
model was chosen for its ability to calculate the interpene- 
trating streams of ions from the conjugate hemispheres of the 
flux tube. 

We realize that there is some discussion as to whether a 
hydrodynamic treatment is valid for plasmaspheric refilling. 
For instance, Schulz and Koons [1972] argued that collision 


frequencies are not nearly high enough to justify a fluid treat- 
ment of even the thermal ions, recommending a kinetic 
approach to modeling plasmaspheric refilling. Performing a 
completely kinetic calculation for the superthermal and ther- 
mal plasma populations is beyond the scope of the present 
study, however, and we must accept the limitations of hydro- 
dynamic modeling for the thermal plasma until a more com- 
prehensive kinetic model is developed. It should be noted, 
however, that Singh et al. [1994] found similar results 
between a hydrodynamic and semikinetic model of plasma- 
spheric refilling for the first hour, and so we will discuss a 
similar timeframe. 

We will discuss the inclusion of the total electron compo- 
nent in the fluid equations (section 2), followed by a descrip- 
tion of the coupling processes in the superthermal electron 
model (sections 3 and 4). Results are presented for the case of 
plasmaspheric refilling along an L = 4 flux tube after a density 
depletion (section 5), closing with a discussion of the results, 
their implications, and further use (section 6). 

2. Total Electron Component in the 
Hydrodynamic Equations 

The hydrodynamic thermal plasma model is described by 
Guiter et al. [1995]. It is a five-moment fluid description for 
the H + and 0 + ions, treating the ions created in the conjugate 
hemispheres as distinct populations (two-stream thermal 
model). The electron density is calculated by assuming quasi- 
neutrality, the electron bulk speed by imposing a zero current 
condition, and the electron temperature from the energy equa- 
tion. These equations are solved along a geomagnetic field 
line under nonsteady conditions. Important modifications 
have been made to this model, however, which will now be 
discussed. 

The basic approach to the thermal electron equations used 
here is the same as before, using the four equations of quasi- 
neutrality, zero current, electron momentum, and electron 
energy to find the four variables n e , u e , E\\, and T e . However, 
the form of these equations used by Guiter et al. [1995] 
assumes that the only electron species is the thermal electron 
population. In a calculation where the electrons are split into 
thermal and superthermal populations, both of these compo- 
nents must be taken into account in all four of these equations. 
In a relatively dense plasma of a filled flux tube, the 
superthermal contribution will be small, but during transient 
events where the thermal plasma population is depleted, the 
superthermal electrons could significantly affect the thermal 
plasma parameters. 

By including superthermal electrons, the condition of quasi- 
neutrality now has the form 

n e +«j =X”' 0) 

f 

and the zero current assumption becomes 

n e u e +n s u s = ^n.u, (2) 

i 

where n s and n s u s are the superthermal electron density and 
flux, respectively, found by taking the appropriate velocity- 
space integral of the superthermal electron flux, <p. 
Examination of (1) and (2) shows that the presence of 
superthermal electron terms will decrease the corresponding 
thermal electron parameters. 
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In deriving the electron moment equations, the electron dis- 
tribution function, /, can be separated into two parts: the 
thermal electrons, f e \ and the superthermal electrons, f s (for a 
complete discussion of this, see Khazanov [1979]). The 
momentum and energy equations then become 

d , \ 2\ 

m e —{ n e u € + n s u s ) + m e B— — \n e u e +n s u s ] 
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where the "es" collision terms represent the momentum and 
energy transfer terms due to Coulomb collisions from the 
thermal electrons to the superthermal electrons. The final 
term in these two equations represent superthermal electron 
losses as particles cascade down to low energies and become 
part of the thermal electron population. Subtracting these 
equations from (3) and (4) yields thermal electron equations 
similar to those of Guiter and Gombosi [1990], 
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where B is the geomagnetic field strength; P a and Q a are the 
pressure and heat flux, respectively, of component a; and the 
collision terms ( 8M e /8t , 8M s /8t f 8E e /8t , and 8E s /8t) rep- 
resent interactions between the electrons and all other species 
(ions and neutral particles). The superthermal electron source 
S s is included in the kinetic equation as a differential source 
term which includes primary, secondary, and tertiary super- 
thermal electron production calculated from an EUV solar spec- 
trum and photoabsorption and impact ionization cross sec- 
tions of the atmospheric species [Khazanov and Liemohn , 
1995]. For conservation of mass, the local production of 
superthermal electrons must balance the local production of 
ions, 

I 

where 5; contributes to the ion moment equations. 

The superthermal electron distribution function, however, 
is calculated by numerically solving the kinetic equation (this 
will be discussed later). We can therefore Find the momentum 
and energy conservation for the superthermal electron compo- 
nent, 


... d(n s u s ) 
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where the last terms in (5) and (6) account for momentum and 
energy transfer from the superthermal electrons to the thermal 
electrons through Coulomb interactions, and have the form 



where E is energy; <p is the superthermal electron flux, 
<p=2Ef/m 2 \ /?=1.7xl0‘ 8 eV 1/2 s cm’ 1 ; A=27te 4 /nA=2.6xl0 12 
eV 2 cm 2 ; and A<p n is the net directional flux along the Field 
line, 

A<p\ \ = 2 njfKpdn 

-1 

where fX is the cosine of the pitch angle, defined as the angle 
between the geomagnetic Field and the particle velocity. The 
energy E min is the lowest energy of the superthermal electron 
distribution where the thermal and superthermal distribution 
functions intersect. (Note that E ( | is the field-aligned electric 
field, whereas E or E a is kinetic energy.) Equation (8) is the 
definition of the thermal electron heating rate due to Coulomb 
collisions with superthermal electrons. 

Using (2) and the thermal ion momentum equations, the 
electron momentum equation (5) can be rewritten and solved 
for En, 
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where p e -m e n e , p e -n e k B T e , and is found from the electron 
energy equation (6), rewritten as a diffusion equation, 
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where Q e has been replaced by the thermal conductivity flux 
according to Banks and Kockarts [1973]. Equation (9) is 
equivalent to the multiple-electron-species electric field 
derived in Mitchell et ai [1990]. Equations (1) and (2) can be 
used to find the thermal electron density and velocity in (9) 
and (10), and these four equations provide the basis for the 
coupling of the two components. 


3. Superthermal Electron Kinetic Equation 
in the Presence of an Electric Field 

A detailed description of the superthermal electron transport 
model is given by Khazanov et ai [1993, 1994] and Khazanov 
and Liemohn [1995], numerically solving the time-dependent, 
guiding center kinetic equation along a geomagnetic field line. 
These studies, however, omitted electric fields. The guiding 
center kinetic equation in the presence of a magnetic field- 
aligned electric field is [cf. Khazanov et ai, 1994, Equation 
(2)] 

P d* * 1 -H*( F , 1 aoa? 

4e dt P ds 2 { E B dsjdfi 

where F represents field-aligned forces such as that from an 
electric field, F--eE\ s is distance along the magnetic field; q 
is the source term; and (S) represents the collision operators 
[e.g., Khazanov and Liemohn ; 1995, Equations (2)-(6)]. 

The electric field alters the magnitude of the parallel veloc- 
ity for a given particle. As can be seen from (1 1), inclusion of 
£|l introduces a drift of the particles in velocity space (drifts in 
energy and pitch angle) as it changes vj|. For instance, if E\\ is 
directed upward, then upward flowing electrons will be deceler- 
ated, shifting to larger pitch angles and lower energies, while 
the downward flowing electrons will be accelerated, shifting to 
smaller pitch angles and higher energies, as shown in Figure 
2. Since both E and fi depend on Vjj, the magnitude of the shift 
in terms of energy and pitch angle is different for each point in 
E-p space. Field-aligned particles will have very little shift in 
pitch angle, while particles with v 1 »vu will have very little 
change in energy. In Figure 2, an electron starting with 
energy E 3 and pitch angle 0j will move to energy E 2 and pitch 
angle 0 2 , while an electron starting at E 2 and 0 } will end up at 



Figure 2. Schematic of the drift in velocity space due to a 
magnetic field-aligned electric field. 


Ej and 0 3 . The energy difference between E } and E 2 is similar 
to that between E 2 and E 3> since both particles started at the 
same pitch angle and the energy dependence is weak, but the 
resulting pitch angles are different because this drift depends 
on both energy and pitch angle. A similar diagram is shown 
for downward flowing electrons. 

Along a geomagnetic field line, the ambipolar electric field 
generally points from each ionosphere toward the equatorial 
plane as plasma is produced in the ionosphere and moves 
upward (Figure 1). Thus electrons are typically decelerated for 
the first half of their plasmaspheric journey, and then acceler- 
ated by roughly the same amount through the second half of 
the plasmasphere. The magnetic field B however, also acts 
upon the particles, focusing them in pitch angle as they move 
toward the minimum field near the equator and then defocusing 
the pitch angle distribution as they move toward the conjugate 
ionosphere. These processes therefore are competing, push- 
ing the electrons in opposite directions in pitch angle. 

Examination of the left-hand side (LHS) of (11) reveals that 
electrons will not move along simple Cartesian grid lines. 
Figure 3 shows typical paths electrons would travel in the s-ji 
plane. They are also drifting up and down in energy, so Figure 
3 is not a level cutaway at a given energy. This drifting as the 
particles move along the field line means that the timescale 
for the distribution function to change in energy and pitch 
angle is less than a plasmaspheric bounce period, r B (the time 
it takes to traverse the field line, mirror in velocity due to the 
inhomogeneous magnetic field, and traverse back to the start- 
ing point). Imposing a Cartesian grid throughout the ( s , E, fi) 
phase space would require a very small time step and a high- 
resolution numerical technique to decrease undesirable compu- 
tational effects associated with approximation errors in the 
d/dE and d/dfi derivatives, making the computation pro- 
hibitively cumbersome. It is desirable, then, to pick a new set 
of variables that would eliminate two of the three drift terms 


LIEMOHN ET AL.: SELF-CONSISTENT PLASMASPHERIC REFILLING 


7527 


Southern latitude base 



Northern latitude base 


Figure 3. Schematic of superthermal electron trajectories in 
space and pitch angle from the LHS of (11), following 
Khazanov et al [1992]. The shaded region is the loss cone, in 
which particles reach the ionospheres, and the striped region 
is the trapped zone, in which particles mirror before the base 
of the flux tube. 


on the LHS of (11), thus having velocity-drift timescales much 
longer than t#. 

To reduce the LHS of (1 1) to spatial advection only, a trans- 
formation of variables must occur, as was performed in the 
previous descriptions of this model (see, for example, 
Khazanov el al. [1993]). The new variables should be chosen 
to eliminate the d/dE and d/dp derivatives on the LHS of (11) 
so that <f> becomes a slowly varying function with s. 
Collisions will then be the only process causing <p to deviate 
from simple advection along the Field line. In general, plas- 
maspheric collisional processes occur on a timescale longer 
than T fi . The transformation presented here includes changing 
the energy and pitch angle variables due to the presence of B 
and a nonzero F. 

We will be transforming from {s, E, p) to ( s , £, p 0 ). The 
new variable set is determined from the characteristics of the 
LHS of (11) by equating the spatial and energy derivative 
terms. 


P -eE\\E\i 


and the spatial and pitch angle derivative terms. 


ds 

V 


dp 

l-// 2 f eEu t 1 dB\ 
2 V E + B ds ) 


(13) 


Using an electrostatic definition of the electric field, F n =- 


d<P/ds , and introducing an electric potential difference, 
A0(s)=0(s)-0(s re f), integration of (12) yields the new vari- 
able £, 


£{s,E) = E-eA<P(s) (14) 

which is the total energy of the particle. Substituting this 
into (13) and integrating reveals the first adiabatic invariant. 


Vo{s, E,m) = ] 1 - 


B n E 




(>v) 


(15) 


where the subscript "o" indicates the parameter value at refer- 
ence point s Q . The potential differences are measured from an 
arbitrary reference point, s re p which does not have to be equal 
to v 

Using (14) and (15), the kinetic equation (11) can be rewrit- 
ten in the form 




where 0'=0U, E, p, s —H, 8, p 0 > s) is the differential flux of 
electrons in the new variable set, Q’ and (5') have also been 
transformed, and E and p now have the form 

E(s,e) = 8 + eA0(s) 


H(s,£,jJ. 0 ) = 



B(s)[£ + eA<P () }i 2 \ 


As seen in (16), the LHS of the kinetic equation is now reduced 
to only advection through the flux tube. 

For this study, the collision terms of Khazanov and 
Liemohn [1995] will be used. While the inelastic collisions 
with neutral particles only require a shift in energy according 
to (14), the Coulomb collision operators and the elastic scat- 
tering with neutrals involve energy and pitch angle transfor- 
mations of the form 


d ( _ d ( J 1^7*0 d(p 

dE\E) d£\E) 2E{£ + eA0o)tlo dp Q 


_d_ 

dpi 


where 




[V ** hn\ 
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n(s)= 


£ + eA&(s) B 0 


£ + eA<P 0 Z?(s) 


and E on the right-hand side is defined above. As before, terms 
of order m/m, and the second derivative with respect to energy 
will be omitted from the kinetic equation calculations 
[Khazanov et al, 1994], and we come to a structure similar to 
that described by Khazanov and Liemohn [1995]. 


4. Numerical Implementation 

It is now possible to define a Cartesian grid throughout ( s , 
£, p a ) phase space and perform the calculation. A typical tra- 
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Figure 4. Schematics of superthermal electron trajectories at constant £ from the LHS of (16). Shaded and 
striped regions have the same meaning as in Figure 3, and white indicates forbidden areas outside of the region 
of existence, where no calculation is necessary. 


jectory in the s-fi 0 plane at constant £ is shown in Figure 4a, 
where the shaded region represents pitch angles that penetrate 
to the ionospheres and the striped region represents the 
trapped zone in which particles mirror and remain in the plas- 
masphere. White shows places outside the region of exis- 
tence, where no calculation will take place. This schematic is 
found from (15), where £ has been chosen larger than 
A<t> 0 +E min . The loss cones (shaded regions) are defined by 
letting s=Sb ase and JJ,= 0 , and is constant for a given £. The 
edge of the trapped zone defining the region of existence is 
found by setting |i-0 for each altitude along the field line. In 
the absence of collisions, the particles travel in straight lines 
in this plane until they reach the edge of the region of exis- 
tence, whereupon they cross over to the stream flowing in the 
opposite direction. This crossover point is the mirror point 
for that particular where fx{fi 0 y s)= 0. The picture is sym- 
metric about the s axis, but not necessarily symmetric about 
the fi 0 axis (this depends on the spatial symmetry of B and F). 
The altitude of the j u Q axis is located at the point of minimum 
7j, defined in (17). As the electrons move from the iono- 
spheres toward the equatorial plane, their kinetic energy 
changes with AO(j) according to (14), even though they 
remain on the same £ level. 

We are particularly interested in simulating the super- 
thermal electrons at low energies, since these particles have 
the largest interaction cross sections with the thermal plasma, 
so it is necessary to choose an £ grid extending down to the 
intersection of the thermal and superthermal electron popu- 
lations. However, in the case of plasmaspheric refilling being 
studied here, A decreases with increasing altitude up to the 
equatorial plane (Figure 2). When leAOCObE-E,™,,, the kinetic 
energy of these particles becomes less than the cutoff E min at 
some point s along the field line, because E=£+eA<&(s). When 
this occurs, the region of existence will end at this altitude for 
this £, and the s-fi 0 plane at constant £ will look like either 
Figure 4b or 4c. In these cases, the reference altitude for ji 0 
has shifted down the field line to the point where tj(s) reaches 
a minimum and E is still greater than or equal to E min . In 


Figure 4b, this point is the altitude when E=E min , while in 
Figure 4c this point occurs lower down the field line, creating 
a trapped population with both mirror points on the same side 
of the equatorial plane. As with Figure 4a, Figures 4b and 4c 
are symmetric about the s axis but not about the )d fl axis. In 
fact, only the s 0 point will have fi 0 defined completely across 
from -1 to +1; for all other altitudes, the region of existence 
will end before reaches zero. 

Figures 4b and 4c represent decoupled hemispheres for 
superthermal electron transport. The particles formed in the 
ionosphere with at an energy that has one of these s-}jl 0 planes 
will not reach the conjugate ionosphere; it will be reflected 
when the potential barrier has removed its field-aligned energy 
and start to move back toward its source ionosphere. This is 
shown in Figures 4b and 4c by the curved gray arrows con- 
necting the upward flowing stream with the downward flowing 
stream at this reflection point. For the purposes of this study, 
the fluxes of particles are simply mapped to the downward 
stream. In Figure 4c the trapped zone does not necessarily 
reach the reflection altitude. This indicates that particles can 
become trapped in one hemisphere, mirroring before the equa- 
torial plane, and usually not much above the ionosphere. This 
trapping mechanism is analogous to the trapped population 
mentioned in the classification scheme of Lemaire and Scherer 
[1972]. That paper discussed the various populations present 
on a polar cap field line, where the ambipolar electric field is 
important in determining ion outflows. Competition between 
the magnetic field divergence and the electric potential barrier 
can sometimes create a trapped electron population with two 
mirror points along an open field line, as seen in the regions 
of existence presented in this study. Chiu and Schulz [1978] 
also discussed this population in the context of parallel elec- 
tric fields along auroral field lines. 

The other aspects of the numerical implementation, initial 
conditions, and boundary values of each model are the same as 
in previous works [Khazanov and Liemohn , 1995; Guiter et 
ai, 1995]. However, the ion heat conduction terms in the 
thermal plasma model are now treated using a fully implicit 
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scheme in order to minimize numerical oscillations which 
occur at shocks when using a Crank-Nicholson scheme. To 
couple the two models, a driver program is used to alternately 
call each model, advancing each solution the same time incre- 
ment and trading the relevant information between the codes. 
The thermal plasma model is started first, advancing half of a 
time increment, followed by a full time increment of the super- 
thermal electron model, so the two models leapfrog in time. 
Using this method, the nonlinearity of the equations due to the 
self-consistent coupling is substituted for the values at the 
previous time step or half time step. The effects of this substi- 
tution are minimized by choosing a small time step. 

5. Results 

The main purpose of this paper is to demonstrate the effects 
of self-consistently collisionally and electrodynamically 
coupling the thermal and superthermal plasmas in a time- 
dependent calculation. This will be shown through plasma- 
spheric refilling after a geomagnetic disturbance along a flux 
tube at L= 4 (geomagnetic equatorial plane crossing at 4 R E 
from the center of the Earth). The simulations start on the 
dayside (0920 LT, 1400 UT) on March 21, 1982. An L= 4 flux 
tube is a closed field line, and so the simulation range extends 
from the base of the northern ionosphere through the plasma- 
sphere to the base of the conjugate ionosphere in the southern 
hemisphere. The conjugate ionospheres are not illuminated 
symmetrically due to the tilt of the dipole, but the source terms 
are within 10% of each other. The initial conditions for the 
thermal plasma along the field line are shown in Figure 1 of 
Gutter et ai [1995], and the calculations begin with no super- 
thermal electron population along the field line. These 
choices are somewhat arbitrary because of a lack of experi- 
mental theoretical knowledge of depleted flux tube character- 
istics, but the use of /t e =0.4 cm' 3 at the equatorial plane is rea- 
sonable [cf. Singh and Horwitz , 1992], 

5.1. Effects of Electrodynamic Coupling 

To demonstrate some of the effects of the electrodynamic 
coupling processes included in this study, a comparison 
between two cases will be shown, with and without the deri- 
vations described above. The results marked "With EDC" 
include the superthermal electron population in the thermal 
plasma equations (section 2), and include the electric field in 
the superthermal electron kinetic equation (section 3). The 
results marked "No EDC" include only the heating rate in the 
electron energy equation but omit superthermal electrons in 
the other thermal plasma equations, and do not include the 
electric field in (1 1). All other aspects of the simulations are 
the same. 

The influence on the protons is demonstrated in Figure 5, 
showing the H + bulk velocity from the southern hemisphere 
stream along the field line after (a) 1 min, (b) 10 min, and (c) 
15 min of refilling. In these plots (as in later figures), dis- 
tance is counted from the base of the northern ionosphere, so 
the northern ionosphere is to the left and the southern iono- 
sphere is to the right, with the equatorial plane located at 
5=4.51 R e . Note that a positive velocity indicates flow from 
north to south, so the large negative velocities mean the 
protons are flowing toward the northern ionosphere. The new 
formulation appears to very slightly accelerate the ions, so 
that by 15 min the shock front is 500 km farther downstream 



Figure 5. Southern stream H + velocities along the field line 
with (solid line) and without (dotted line) electrodynamic 
coupling at (a) 1 min, (b) 10 min, and (c) 15 min of refilling. 
The dashed line shows w=0. Distances are counted from the 
base of the northern ionosphere. 


than without the new coupling terms. The maximum speed of 
the shock is also 1 km/s higher with the new terms at 15 min. 
This additional progress of the ions can be explained by the 
influence of the superthermal electrons on the electric field 
during the initial stages of the refilling process. 

The impact of the superthermal electron terms on the 
thermal electron population is significant. Figure 6 shows the 
thermal electron velocity as a function of field-aligned dis- 
tance at several times, with and without the electrodynamic 
coupling processes. After 10 s of refilling time (Figure 6a), 
the electrodynamic coupling is causing the thermal electrons 
to stream away from the equator, whereas the thermal electrons 
stream toward the equator when this process is not included. 
This is due to the large superthermal electron flux now taken 
into account in the currentless condition (2). A result such as 
this could be unstable, but the topic of stability in the thermal 
and superthermal distributions will be addressed at a later date. 
After 1 min (Figure 6b), the results with electrodynamic 
coupling are beginning to resemble the results without, as the 
ions flow into the plasmasphere and balance the superthermal 
electron flux. By 15 min of refilling, as seen in Figure 6c, the 
two results are quite similar, with the only deviations being 
near the ionospheres. This is because the counterstreaming 
superthermal electron populations are not balanced at these 
altitudes, because one stream is near its source and the other 
stream has traversed the plasmasphere. 
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Figure 6. Thermal electron velocities along the field line 
with (solid lines) and without (dotted lines) electrodynamic 
coupling at (a) 10 sec, (b) 1 min, and (c) 15 min of refilling. 
The dashed lines show u= 0. 


From section 3, the use of (16) instead of (11) leads to 
changes in the super-thermal electron results. Figure 7a shows 
omnidirectional flux spectra versus kinetic energy after 10 s of 
refilling at the topside ionosphere and the equatorial plane, 
with and without electrodynamic coupling. Omnidirectional 
flux is defined as 



-l 


The results at 800 km clearly show the photoelectron pro- 
duction peaks in the 20-30 eV range, and it can be seen that 
there is a downward shift in the energy of these peaks of 
around 1 eV in the results with the coupling terms. At 4.51 
/?£, the potential drop is close to 5 V, and the entire distri- 
bution has shifted downward accordingly. An increase at low 
energies due to the coupling is evident in the equatorial 
results, also due to the downward shift of the distribution 
function. Coulomb collisions, which have an energy depen- 
dence of E ‘ 2 , are more efficient at degrading the superthermal 
flux levels at lower energies, and so the results with electrody- 
namic coupling are not a perfectly shifted image of the results 
without coupling. 

Figure 7b shows the effects of electrodynamic coupling on 
the superthermal electron pitch angle distributions. Presented 
here are equatorial distributions after 30 min with and without 
coupling at 5 and 10 eV of kinetic energy. The fluxes that 
include electrodynamic coupling processes are lower than the 
results without at small pitch angles, but are higher at larger 


pitch angles. This effect is due to the decrease in vj| from the 
field-aligned electric field, causing a shift to lower energies 
and to larger pitch angles as the electrons move from the 
ionosphere to the equatorial plane (see Figure 2). The effect is 
small because the superthermal electron source is the iono- 
sphere. In the plasmasphere, then, the source is in the loss 
cone, located at small pitch angles, where the vjj influence is 
primarily in energy and not in pitch angle (see equation (11) 
and Figure 2). Coulomb collisions with the thermal plasma 
are responsible for scattering electrons into the trapped zone, 
and the electric field acts to push them further into the trapped 
zone. For reference, the conventional loss cone (5 influence 
only) for an L= 4 flux tube is 5.5°. Here, however, the nonzero 
F term in (11) expands the loss cone as a function of energy 
and time. However, even for a case such as £= 6 eV and AO 0 =-5 
V, the loss cone will only expand to 13.5°. 

5.2. Time Dependence of the Results 

It is interesting to show the development of the solution 
with respect to time. As mentioned in the Introduction, this is 
the first model to self-consistently couple the thermal plasma 
and superthermal electron calculations collisionally and elec- 
trodynamically under nonsteady conditions, and here we will 
demonstrate the time dependence of the results. 

With the inclusion of electrodynamic coupling, the total 
ion density must now equal the thermal and superthermal 
electron densities combined. In this simulation, the initial 



(a) 0omni (10 sec) 



Figure 7. Superthermal electron (a) omnidirectional flux 
spectra and (b) equatorial pitch angle distributions with (solid 
lines) and without (dotted lines) electrodynamic coupling. The 
$omni are shown after 10 s of refilling at 800 km (topside 
ionosphere) and 4.51 R E (equatorial plane) along the field 
line. Pitch angle distributions are shown after 30 min of 
refilling at kinetic energies of 5 and 10 eV. 






LIEMOHN ET AL.: SELF-CONSISTENT PLASMASPHERIC REFILLING 


7531 



0 5 10 15 20 25 30 

Time (min) 

Figure 8. Refilling time dependence of the equatorial values 
of (a) electron densities, (b) thermal electron temperatures, and 
(c) electrostatic potential differences with (solid lines) and 
without (dotted lines) electrodynamic coupling. 


thermal electron density at the equatorial plane is 0.4 cm' 3 . It 
is conceivable, then, for the superthermal electron density to 
approach this level and significantly affect the thermal 
electron density. In Figure 8a, equatorial densities for the 
thermal and superthermal components of the electron popu- 
lation are shown for the first 30 min of refilling. For the 
results with the new coupling terms, the superthermal density 
has risen to 0.18 cm' 3 and the thermal density has dropped to 
0.26 cm' 3 after 2 min of refilling. This illustrates the impor- 
tance of self-consistently including superthermal electrons 
during events where the thermal plasma is depleted. After 30 
min, however, the thermal plasma is 2 orders of magnitude 
larger than the superthermal density, and the influence of 
superthermal electrons on the quasi-neutrality condition is 
greatly diminished. The superthermal densities without 
coupling are higher than those with coupling; this is because 
of deceleration due to the electric field. Also notice that the 
superthermal electron density is dropping slightly after 15 
min. This is due to the increased thermal plasma density in the 
plasmasphere increasing the Coulomb losses of the super- 
thermal electron population. 

Figure 8b shows the thermal electron temperature at the 
equatorial plane. The use of the Spitzer conductivity for the 
thermal heat flux is responsible for a rapid decline in the 
plasmaspheric temperature. The coupling results have a lower 
temperature because of the convection shown in Figure 6 
transporting additional energy out of the plasmasphere. The 
temperatures start to increase, very slightly, after 20 min or so 


of refilling. This is consistent with the decrease in super- 
thermal electron density and is because of the Coulomb col- 
lision heating of the thermal electrons, defined in (8). 

In Figure 8c, AO is shown for the first 30 min of refilling at 
the equatorial plane. Here, the reference altitude for the poten- 
tial is taken at the point of maximum O, located in the north- 
ern ionosphere near the F 2 peak. Notice that A0 0 changes 
from -5.2 to -2.2 V during the first 15 min. The initial Ad> 0 
without the superthermal electron influence is -4.9 V, so the 
superthermal electrons increase the potential drop by roughly 
5% during the initial stages of the refilling process. The 
crossover of these two curves is due to the difference in elec- 
tron temperature shown in Figure 8b. Initially, the coupling 
increases the potential drop, but as the thermal electrons cool, 
the coupling has the opposite effect. 

Plate 1 shows the development of the superthermal electron 
distribution function along the field line by presenting omni- 
directional fluxes after 10 s, 2 min, and 30 min of refilling. In 
Plate la, the photoelectron production peaks between 20 and 
30 eV are clearly discernible, and they reveal the downward 
shift of the particles in energy as they move toward the equato- 
rial plane. Plate lb also has the production peaks, but they are 
beginning to be smoothed away by the increasing super- 
thermal electron flux level around them, by Coulomb losses as 
the thermal density refills from the ionosphere, and by the 
decreasing the potential difference. These peaks are barely 
noticeable after 30 min in Plate lc, where the increased 
Coulomb losses act to smooth the peaks out of the distri- 
bution. 

Plate la appears to be a butterfly distribution, except this is 
the E-s plane rather than v„-v x . There are several reasons for 
the formation of this distribution. At low energies, the elec- 
trons have not yet reached the equatorial region. A 3 eV 
electron moves at 1000 km/s, and even if its pitch angle is 
zero, it will penetrate less than 2 Re along the field line into 
the plasmasphere after 10 s. However, Coulomb collisions are 
occurring which help to decelerate the faster particles while 
they are in the plasmasphere, and so the equatorial region is 
not totally devoid of low-energy superthermal electrons. The 
production peaks in the 20-30 eV range are formed by the 
strong He II-30.4 nm solar emission line, which has an energy 
of 40 eV. The ionization potentials of the atmospheric neutral 
particles create spikes in the electron distribution, and these 
electrons can escape the ionosphere and travel through the 
plasmasphere. At higher energies, the streams of super- 
thermal electrons from the two ionospheric connection points 
of the flux tube have interpenetrated and started flowing down 
toward the conjugate ionosphere. However, Coulomb col- 
lisions are very slow to scatter these electrons into the trapped 
zone, and the electrodynamic influence is also very weak. 
Therefore the trapped zone is relatively empty for these ener- 
gies, and since omnidirectional flux is an average over pitch 
angle, the flux appears to decrease in the equatorial region. 
Thus an E-s plane butterfly distribution arises. 

Let us now examine the superthermal electron equatorial 
pitch angle distribution development. Plate 2 shows the time 
dependence of the 5 and 30 eV kinetic energy electrons. Note 
that the two plots have different color scales and dynamic 
ranges. The loss cone fluxes appear to achieve a steady state 
level within a minute since these pitch angles are connected to 
the source regions. The trapped zone, however, must have par- 
ticles scattered into it from the loss cone, and it is clear that 
this process takes time to complete. It is evident that the 5 eV 
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Plate 1* Superthermal electron omnidirectional flux spectra along the field line at (a) 10 s, (b) 2 min, and (c) 
30 min of refilling time. Distances are measured from the base of the northern ionosphere. All three plots 
share the color scale to the right. 
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Plate 2. Superthermal electron equatorial pitch angle distributions vs. refilling time at kinetic energies of 
(a) 5 eV and (b) 30 eV. Note that each plot has its own color scale. 
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Plate 3. Superthermal electron omnidirectional flux spectra along the field line for nonsymmetric illumi- 
nation at (a) 10 s, (b) 2 min, and (c) 30 min of refilling time. Distances are measured from the base of the 
northern ionosphere. Both plots share the color scale to the right. 
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Figure 9. Superthermal electron (a) equatorial pitch angle 
distributions and (b) omnidirectional flux spectra after 10 s 
(dashed lines) and 15 min (dotted lines) of refilling time and 
steady state (solid lines) results. Pitch angle distributions are 
at 10 eV, and 0 omm are at 800 km (topside ionosphere) and 
4.51 Re (equatorial plane) up the field line. 


fluxes are refilling much faster than the 30 eV fluxes because 
of the E' 2 energy dependence of Coulomb collisions. The 
pitch angles near the loss cone are reaching a steady level, but 
notice that the distributions deeper in the trapped zone are still 
refilling after 30 min, and with the thermal plasma still on the 
rise at this point, electrons will continue to scatter to larger 
pitch angles for many more hours. 

5.3. Steady State Results 

The process of refilling and the long time period required to 
reach a steady state level can be accelerated by conducting a 
simulation that starts with a flux tube filled with thermal 
plasma. In this case, only the superthermal electron distri- 
bution needs to develop, and Khazanov and Liemohn [1995] 
discussed this jump to steady state for the superthermal 
electron model. The time step is set to infinity and the model 
iterates to a converged solution. The two models are still elec- 
trodynamically coupled, but the time steps are now iteration 
steps. The thermal plasma density at the equator is taken to be 
2500 cm -3 , representing a filled plasmaspheric flux tube. All 
other aspects of the calculation are the same. 

Figure 9a shows 10 eV equatorial pitch angle distributions 
after 10 s and 15 min of refilling and at steady state. Notice 
that the steady state results have a lower flux at small pitch 
angles but a higher flux at large pitch angles. This isotro- 
pization of the distribution is due to the increased Coulomb 
scattering cross section in the steady state results. Comparing 
this steady state distribution with the 5 and 30 eV distri- 


butions developing in Plate 2, it can be seen that the trapped 
zone will take a long time to reach this steady-state distri- 
bution. 

Steady-state omnidirectional fluxes are also compared with 
the refilling results, shown in Figure 9b. This is similar to 
Figure 7a, showing spectra at 800 km and 4.51 /?£up the field 
line. Notice that the production peaks are smoothed away in 
the steady state results by interactions with the thermal 
plasma. Coulomb collisions have also eroded the low-energy 
end of the steady state results, and in general decreased the flux 
levels compared to the results after 15 min of refilling. The 
potential difference from ionosphere to equatorial plane is 
-1.66 V for the steady state case, so there is less than a volt 
difference between this and the potential difference at 15 min. 

5.4. Source Term Effects 

So far, all of the results have been "symmetric," that is, 
both ionospheric footprints have been illuminated and 
refilling has occurred from both hemispheres in roughly equal 
proportions. It is useful to examine the case with nonsymmet- 
ric illumination: one sunlit ionosphere and the other in dark- 
ness. This will demonstrate the influence of an ionosphere on 
its conjugate point. 

In Plate 3, omnidirectional fluxes along the field line after 
10 s, 2 min, and 30 min of refilling are shown without illumi- 
nation in the southern hemisphere. Compare these results to 
Plate 1 (with southern illumination). Plate 3a shows a dra- 
matic decrease in the flux of low-energy electrons in the sec- 
ond half of the plasmasphere. However, particles are still 
traveling to the conjugate hemisphere, and can deposit their 
energy to the thermal plasma there. Plate 3b is beginning to 
resemble Plate lb, and Plate 3c is very close to Plate 3c. It is 
still apparent, though, that the fluxes in the southern 
hemisphere are lower, but it is also clear that plasma from the 
northern hemisphere flows down to the southern hemisphere 
to interact with the nonsunlit atmosphere. 

Figure 10 shows thermal plasma heating rates along the 
field line. The refilling results, both with and without south- 
ern hemisphere illumination, are shown after 15 min of 
refilling. Also shown are heating rates from the steady state 



state case (solid line) and after 15 min of refilling with (dotted 
line) and without (dashed line) southern hemisphere illumi- 
nation. 
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calculation described in section 5.3. Comparing these results 
to the heating rates shown by Liemohn and Khazanov [1995] 
reveals that the steady state calculation is quite similar to the 
previous steady state results for L= 4. Notice that the non- 
symmetric results are less than the symmetric results for most 
of the field line, dropping orders of magnitude less in the 
southern ionosphere. This is consistent with the depleted 
ionosphere at that end of the flux tube, since after 15 min the 
H + stream from the northern hemisphere is just reaching the 
southern ionosphere. 

Heating of the thermal plasma by the superthermal 
electrons is a slow process during refilling, as can be seen by 
the factor of 30 difference in the heating rate at the equatorial 
plane. This is due not only to the increased thermal plasma 
density but also from the increased superthermal flux in the 
trapped zone. 

5*5. Data Comparison 

Finally, it is beneficial to show a comparison of results 
from this model with experimental data. Figure 1 1 shows 
omnidirectional fluxes from Atmospheric Explorer E (AE-E) 
and steady state model results for similar conditions. The data 
are reproduced from Doering et al. [1976], for day 355 of 1975 
at 182 and 365 km altitude. The solar zenith angles for the 
two spectra at 50° and 37°, respectively. Since AE-E flew in a 
nearly equatorial orbit, the model comparisons are made at 0° 
geographic latitude, choosing a local time appropriate with 
the given solar zenith and assuming the data collection 
occurred in the morning. Initial profiles for the thermal 
plasma are taken from the IRI model [Bilitza, 1990] for these 
conditions in order to start with realistic equatorial iono- 
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Figure 11. Comparison of model results (solid lines) with 
AE-E data (dotted lines) at (a) 182 km and (b) 365 km altitude 
on day 355 of 1975. The satellite data are reproduced from 
Doering et al. [1976]. 


spheric densities and temperatures. In Figure 1 1 a, the spectra 
agree very closely for most of the energy range. The model 
predicts a slightly higher flux in the 5-15 eV range, but this 
difference is less than a factor of two. Figure lib also shows 
good agreement, with the model predicting more definition in 
the 20-30 eV range and lower fluxes above 30 eV by a factor of 
less than two. These differences could be explained by uncer- 
tainties in the experimental data, differences in the neutral 
atmosphere or ionospheric plasma profiles, or uncertainties in 
the collisional cross sections used in the model. The larger 
fluxes at low energy and the increased definition of the pro- 
duction peaks in the model results indicates that the thermal 
plasma density from IRI is lower than the actual densities; a 
higher plasma density would act to smooth out these features 
of the distribution function. The comparison does show, how- 
ever, that the model accurately calculates the main features of 
the photoelectron spectrum. 

6. Discussion 

In this paper, we have discussed the need for simultaneous 
model calculations of the thermal plasma and superthermal 
electrons through self-consistently coupling both collisional 
and electrodynamic processes between the populations. This 
was achieved by introducing the total electron population 
(thermal plus superthermal) into the quasineutrality condition, 
currentless assumption, electron momentum equation, and 
electron energy equation. Knowing the superthermal electron 
distribution from numerically solving the field-aligned, 
guiding center kinetic equation, the former equations are used 
to find the thermal electron density, bulk velocity, and tem- 
perature, as well as the parallel electric field. These quantities 
are then used in the solution of the kinetic equation to obtain a 
new superthermal electron distribution. Of special interest is 
that the ambipolar electric field was included in the kinetic 
equation, necessitating a new variable transformation to 
reduce numerical difficulties arising from velocity-space drifts 
due to the electric and magnetic fields. 

Results from calculations based on this new formalism were 
described in detail. It was shown that the potential drop 
decreases from -5 to -2 V once the initial ion shock fronts 
penetrate to the conjugate ionosphere. This drop in the poten- 
tial is due to the increase in thermal plasma density, decrease 
in thermal plasma temperature, and decrease in superthermal 
electron flux. From there begins the slow refilling of the 
plasmasphere, which will take a minimum of several more 
hours. Removing ion production from one of the ionosphere 
shows that the sunlit ionosphere can still deposit a significant 
amount of energy into the dark hemisphere. 

In the works of Khazanov et al. [1993] and Liemohn and 
Khazanov [1995], it was shown that a depleted flux tube could 
take several hours to reach a steady state level. Here we con- 
firm these results but suggest that superthermal electron 
refilling will progress faster due to the effects of the field- 
aligned potential. These effects include decelerating the 
electrons as they move towards the equatorial plane and push- 
ing them to larger pitch angles. The pitch angle drift directly 
contributes to enhanced plasmaspheric trapping, and the 
energy drift indirectly contributes by moving the particles to a 
lower energy where Coulomb collisions with the thermal 
plasma will have a greater scattering cross section and thus the 
capability of trapping more particles. These trends were 
pointed out in Figure 7, where the coupled results showed 
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enhanced trapping of the superthermal electrons. This effect 
is small, however, and refilling will still take several hours to 
complete. The decreased superthermal electron density at the 
equatorial plane seen in Figure 8a is due to the decreased loss 
cone flux from the potential drop. As the potential drop 
relaxes to smaller values, these loss cone fluxes will return 
since they are constantly replenished by the ionospheric 
source. 

Influences on the thermal plasma were also calculated. The 
results indicate that the thermal plasma will refill faster than 
before due to slightly enhanced proton velocities in the initial 
shock front. This effect is small because of the large inertia 
involved in changing the ion distributions. The dramatic 
influence on the thermal electron velocities shown in Figure 
6a has almost disappeared by Figure 6c. This indicates that 
the thermal plasma will be influenced by the superthermal pro- 
cesses during the early stages of refilling, causing an initial 
enhancement in the ion population and an initial depletion of 
the the thermal electron population, but eventually the 
thermal electron population recovers and refilling will con- 
tinue as before with the thermal population slightly increased 
over previous results due to the initial boost. More substan- 
tial changes would be expected along auroral or polar cap field 
lines where the superthermal electron fluxes would not be bal- 
anced by particles from the conjugate ionosphere, allowing a 
longer timeframe for influences on the thermal ions to 
develop. 

As mentioned in the Introduction, we realize a hydrody- 
namic description of the thermal plasma may not be entirely 
valid for the low densities involved with the early stages of 
plasmaspheric refilling. The initial flow of ionospheric 
plasma through the flux tube will only fill the loss cones of 
the velocity distribution, resulting in counterstreaming flows. 
This situation could be unstable, and Schulz and Koons [1972] 
argued that this two-stream instability is more likely to refill 
the high pitch angle trapped zone of velocity space than either 
Coulomb collisions or collisionless shocks. Indeed, they 
argue against the appearance of a collisionless shock. It 
would be ideal to have a fully kinetic model of the thermal and 
superthermal populations, making no distinction between 
them and modeling the small-angle scattering and possible 
wave generation and interaction. This model is a step in this 
direction by including the ambipolar electric field in the 
kinetic calculation, and we plan to extend the kinetic portion 
of the model down into the thermal energy regime. 
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